      SUBROUTINE SOLMDT(TSEC50,SUNVEC,IFLAG)
C
C     CALCULATE THE POSITION OF THE SUN IN MEAN OF DATE COORDINATES.
C
C
C     VARIABLE    TYPE    I/O         DESCRIPTION
C     --------    ----    ---         -----------
C
C     TSEC50       R*8     I     INPUT TIME IN SECONDS SINCE JAN 1 1950
C                                0.0 HOURS ET. MAY BE ANY VALUE + OR -.
C                                ACCURACY IS NOT HIGH ENOUGH TO WARRANT
C                                DOING AN ET-UT CORRECTION. TIME CAN BE
C                                UT IF IT IS MORE CONVENIENT.
C
C     SUNVEC(3)    R*8     O     VECTOR TO THE SUN FOR THAT TIME
C
C     IFLAG        I*4     I     EQ 1, SUNVEC IS A UNIT VECTOR
C                                NE 1, SUNVEC IS IN KM.
C
C***********************************************************************
C
C  BY C PETRUZZO. ABOUT 1/83. COMMENT AND MINOR CODE MODS 7/83.
C
C***********************************************************************
C
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION SUNVEC(3)
      DATA TWOPI/6.283185307179856D0/,SECDAY/86400.D0/
      REAL*8 DEGRAD/ 57.29577951308232D0 /
C
C  SEE EXPLANATORY SUPPLEMENT TO THE ASTRONOMICAL EPHEMERIS AND THE 
C  AMERICAN EPHEMERIS AND NAUTICAL ALMANAC(1974,=1961+AMMENDMENTS),
C  PAGE 98. GSFC LIBRARY NUMBER IS QB8.U5, E96.
C
C  ACCURACY: THE RESULTS OF USING THIS ROUTINE WERE COMPARED AGAINST
C            THE 1983 ASTRONOMICAL ALMANAC TABLE VALUES FOR RIGHT
C            ASCENSION AND DECLINATION. COMPARISONS WERE DONE AT
C            0.0 HRS ON 1/1/83 AND AT 10 DAY INTERVALS THROUGHOUT 1983.
C            THE MAXIMUM ERROR FOUND WAS ABOUT 2 ARC MINUTES.
C
      INTEGER INIT/1/
C
C  INITIALIZE PARAMETERS TO VALUES AT JAN 0.5, 1900, E.T.
      IF(INIT.EQ.1) THEN
        INIT=0
        GAM0 = 281.220844D0/DEGRAD
        GAMD = .0000470684D0/DEGRAD
        G0 = 358.475833D0/DEGRAD
        GD = 0.985600267D0/DEGRAD
        EPS0 = 23.452294D0/DEGRAD
        EPSD = 0.35626D-06/DEGRAD
        END IF
C
C
C  CONVERT FROM SECONDS SINCE 1/1/50, 0.0 HRS E.T., TO DAYS SINCE 
C  JAN 0.5, 1900, E.T.
      D = TSEC50/SECDAY + 18262.5D0
C
C  ECCENTRICITY OF EARTH'S ORBIT ABOUT THE SUN AT THE INPUT TIME.
      E = 0.1675104D-01 - 0.11444D-08*D
C  OBLIQUITY OF THE ECLIPTIC AT THE INPUT TIME.
      EPS = EPS0 - D*EPSD
C  ANGLE FROM EQUINOX TO PERIHELION AT THE INPUT TIME.
      GAM = GAM0 + D*GAMD
C  MEAN POSITION OF THE SUN, REFERRED TO PERIHELION, AT THE INPUT TIME.
      G = G0 + D*GD
C  LONGITUDE OF THE SUN AT THE INPUT TIME.
      AL = GAM + G + 2.D0*E*DSIN(G)
      AL = DMOD(AL,TWOPI)
C
C  CONVERT ANGLES TO UNIT VECTOR.
      SINAL = DSIN(AL)
      SUNVEC(1) = DCOS(AL)
      SUNVEC(2) = SINAL*DCOS(EPS)
      SUNVEC(3) = SINAL*DSIN(EPS)
C
C  IF WANTED, CONVERT UNIT VECTOR TO POSITION VECTOR.
      IF(IFLAG.NE.1) THEN 
        R=149597871.D0/(1.D0+E*DCOS(G))
        SUNVEC(1)=R*SUNVEC(1)
        SUNVEC(2)=R*SUNVEC(2)
        SUNVEC(3)=R*SUNVEC(3)
        END IF
C
      RETURN
      END
